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ABSTRACT 

In  this  thesis,  an  adaptive  two  dimensional  least  mean  squares  (LMS)  algorithm  and 
a  recursive  least  squares  (RLS)  algorithm  are  developed  from  the  one  dimensional  algo- 
rithms. Design  of  the  two  dimensional  LMS  and  RLS  algorithms  are  studied  for  accu- 
racy based  on  the  results  of  a  two  dimensional  system  identification  model  which  was 
used  in  testing  the  algorithms.  Application  of  the  two  algorithms  is  demonstrated 
through  computer  simulation  in  which  the  adaptive  filters  are  employed  in  a  noise 
canceler  and  an  adaptive  line  enhancer  and  applied  to  an  image  processing  problem. 
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I.     INTRODUCTION 

The  area  of  digital  signal  processing  has  experienced  a  rapid  growth  in  the  last  dec- 
ade. Reasons  for  this  have  been  the  tremendous  advances  in  integrated  circuit  technol- 
ogy, and  some  significant  developments  in  digital  processing  techniques  achieved  during 
this  period.  Included  in  these  developments  are  methods  v^^'hich  extend  certain  one  di- 
mensional digital  signal  processing  techniques  to  two  dimensions.  This  extension  is  by 
no  means  trivial.   Three  significant  factors  must  be  considered: 

1.  more  degrees  of  freedom  are  available  in  a  two  dimensional  system;  this  gives  a 
system  designer  more  flexibility  than  the  one  dimensional  case, 

2.  one  dimensional  problems  generally  involve  considerably  less  data  than  two  di- 
mensional problems,  and 

3.  the  mathematical  methods  for  handling  two  dimensional  systems  are  generally  less 
complete  than  those  for  one  dimension. 

As  the  techniques  for  processing  multidimensional  data  have  improved,  the  appli- 
cations of  digital  signal  processing  have  spread  from  one  dimensional  to  two  dimen- 
sional, to  n-dimensional  data.  There  exists  many  physical  phenomena  that  inherently 
depend  on  two  or  more  independent  variables.  In  the  prediction  of  weather  and  in 
seismic  analysis,  the  data  depends  on  more  than  one  independent  variable.  Moreover, 
data  originating  from  one  dimensional  processes  can,  in  some  cases,  be  considered  two 
dimensional.  For  instance,  data  from  periodic  or  cyclic  processes  can  be  represented  as 
two  dimensional  arrays  by  using  their  periodicity. 

Besides  the  general  applicability  of  two  dimensional  signal  processing  in  the  above 
cases,  other  areas  which  have  experienced  significant  growth  in  recent  years  include  ra- 
dar, sonar,  and  radio  astronomy.  The  two  dimensional  processing  of  images  is  also  a 
very  important  one.  Images  depend  on  two  spatial  variables,  and  are  continuous  in 
nature.  However,  if  we  digitize  them  and  assume  linear  models  in  their  formation,  dis- 
tortion, and  recording  we  then  have  techniques  that  can  be  used  in  their  processing. 
Depending  on  the  applications,  diflerent  processing  techniques  are  used,  notably:  en- 
hancement and  restoration  of  images,  and  segmentation  and  encoding  of  images.  For 
details,  see  References  1,2,3,  4  . 

This  brief  discussion  about  applications  of  two  dimensional  digital  signal  processing 
shows  that  they  are  found  in  a  wide  variety  of  fields.    In  this  thesis,  we  are  interested  in 


extending  one  dimensional  adaptive  filtering  techniques  to  two  dimensions  and  then 
applying  this  in  the  area  of  image  restoration.  Once  the  transition  from  one  to  two  di- 
mensions is  understood,  the  extension  to  n-dimensional  signal  processing  is  fairly 
straightforward. 

A.  OBJECTIVES  OF  THE  THESIS 

The  use  of  the  Wiener  filter  has  proven  to  be  a  very  powerful  tool  in  the  area  of 
image  restoration  [Refs.  1,2,  4  ].  However,  one  disadvantage  is  the  fact  that  the  Wiener 
filter  operates  under  the  assumption  that  the  image  is  stationary  which  generally  is  not 
the  case  in  an  image  processing  problem.  In  one  dimensional  signal  processing  the 
Wiener  filter  concept  has  provided  a  basis  for  various  adaptive  filtering  algorithms. 
Within  these  adaptive  algorithms,  the  filter  possesses  characteristics  which  can  be  mod- 
ified to  achieve  some  end  or  objective  and  is  usually  assumed  to  accompHsh  this  "adap- 
tation" automatically  without  the  need  for  substantial  intervention  by  the  user.  The 
adaptive  filter  can  "learn"  the  signal  characteristics  when  first  turned  on  and  thereafter 
can  track  changes  in  these  characteristics.  The  first  objective  of  this  thesis  is  to  examine 
a  two  dimensional  Wiener  model  for  image  restoration  which  can  then  be  extended  to 
adaptive  algorithms.  Due  to  its  relatively  low  computational  requirements  and  the  fact 
that  it  will  work  in  a  variety  of  signal  environments,  the  least  mean  square  (LMS)  algo- 
rithm will  be  investigated  first. 

The  second  objective  of  this  thesis  is  to  develop  a  second  two  dimensional  adaptive 
algorithm.  In  this  case,  we  will  work  with  the  recursive  least  square  (RLS)  algorithm 
which  offers  faster  convergence  than  the  gradient-search-type  algorithms  but  generally 
involves  a  greater  cost  per  data  sample  and  more  numerical  difficulties. 

The  final  objective  is  to  implement  the  LMS  and  RLS  algorithms  within  a  system 
identification  model,  a  noise  canceler,  and  an  adaptive  line  enhancer.  Each  algorithm 
possesses  several  variants  and  various  parameters  which  can  be  modified.  A  represen- 
tative sample  of  the  various  outputs  will  be  examined  and  the  results  compared  in  order 
to  see  which  provides  the  optimum  solution  under  given  conditions. 

B.  ORGANIZATION  OF  THE  THESIS 

Chapter  II  is  designed  to  review  the  development  of  the  one  dimensional  Wiener 
filter  and  then  extend  it  to  two  dimensions  where  it  can  be  incorporated  into  a  two  di- 
mensional LMS  adaptive  algorithm.  Computer  simulation  of  the  algorithm  within  a 
system  identification  model  is  performed  and  the  results  are  shown.  The  two  dimen- 
sional RLS  algorithm  is  derived  in  Chapter  III  and  again  the  results  of  the  algorithm 


within  a  computer  simulated  system  identification  model  are  shown.  Chapter  IV  con- 
tains the  results  of  implementing  the  LMS  and  RLS  algorithms  in  a  noise  canceler  and 
an  adaptive  line  enhancer.   Conclusions  concerning  the  results  are  also  presented. 


II.     TWO  DIMENSIONAL  ADAPTIVE  LEAST  MEAN  SQUARE 

ALGORITHM 

In  this  chapter,  we  will  review  the  derivation  of  the  one  dimensional  Wiener  filter 
and  then  develop  an  algorithm  to  extend  it  to  the  two  dimensional  case.  We  will  use  the 
two  dimensional  Wiener  filter  to  achieve  a  two  dimensional  least  mean  square  (LMS) 
adaptive  filter  algorithm  which  will  be  demonstrated  to  be  useful  in  image  processing  by 
applying  the  algorithm  in  a  noise  canceler  mode  and  in  an  adaptive  line  enhancer  con- 
figuration for  the  restoration  of  images  corrupted  by  noise. 

A.     ONE  DIMENSIONAL  WIENER  FILTER 

The  problem  of  estimating  one  signal  from  another  is  one  of  the  most  important  in 
signal  processing.  In  many  applications,  the  desired  signal  is  not  available  or  observable 
directly.  Instead,  the  observable  signal  is  a  degraded  or  distorted  version  of  the  original 
signal.  The  signal  estimation  problem  is  to  recover,  in  the  best  way  possible,  the  desired 
signal  from  its  degraded  replica.  One  typical  example  which  we  will  deal  with  in  this 
paper  is  an  image  recorded  by  an  imaging  system  that  has  been  corrupted  by  noise.  The 
problem  is  to  undo  the  noise  induced  distortion  and  restore  the  original  image. 

This  represents  the  classic  one  dimensional  problem  in  communication  theory-  where 
we  must  obtain  an  estimate  of  a  signal  of  interest,  which  can  be  observed  in  the  presence 
of  some  additive  noise.  In  other  words,  the  available  information  about  the  signal, 
s{n),  is  contained  in  the  received  signal: 

u{n)  =  sin)  +  win)  (2.1) 

where  win)  is  the  noise.    Therefore,  we  must  process  this  available  signal  w(«)  through 
an  optimal  processor  that  produces  the  best  possible  estimate  of  sin). 

In  order  to  establish  a  two  dimensional  Wiener  filter,  we  must  first  develop  a  one 
dimensional  algorithm.  This  task  has  been  approached  from  many  different  directions 
[Refs.  5,6,7].  The  formulation  by  Haykin  [Ref  8]  provides  the  most  logical  extension  to 
two  dimensions.  First,  we  consider  a  tapped  delay  line  filter  similar  to  Figure  1  on  page 
5.  The  filter  consists  of  a  set  of  delay  elements,  a  corresponding  set  of  adjustable  tap 
gains  or  coefilcients  /i(l), /2(2),...,  A(i\/)  connected  to  the  tap  inputs,  and  a  set  of  adders 
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Figure  1.      Tapped  Delay  Line  Filter 

for  summing  the  resultant  outputs.   The  filter  is  driven  by  a  random  time  series  produc- 
ing the  sequence  u{n),  u{n  —  1),...,  u{n  —  iV/+  1)  as  the  M  tap  inputs  of  the  filter. 

We  can  express  the  signal  produced  at  the  filter  output,  y{n),  by  the  convolution 
sum: 


y{n)  =  Yh{k)u{n  -  A  +  1). 


(2.2) 


We  desire  a  filter  which  in  some  way  minimizes  the  difierence  between  some  desired  re- 
sponse, d{n),  and  the  corresponding  value  of  the  actual  filter  output.  This  difierence  can 
be  denoted  as 


e{n)  =  d{n)  -  y{n) 


(2.3) 


where  e{n)  is  called  the  error  signal.    In  Wiener  theor}',  the  filter  is  optimized  by  mini- 
mizing the  mean-square  value  of  this  error  signal,  e{n)  . 


Let  the  mean-square  value  of  the  error  be  denoted  by 

MSE=E{e\n)}  (2.4) 

where  E  {.}  is  the  expectation  operator.  This  mean-square  value  is  a  real  and  positive 
scalar  quantity,  representing  the  average  normalized  power  of  the  error  signal,  e{n). 
Substituting  Equation  (2.3)  into  (2.4)  yields 

MSE  =  E  {d\n)}  -  2  £  {d{n)y{n)}  +  E  {y\n)}.  (2.5) 

Next,  substituting  Equation  (2.2)  into  Equation  (2.5)  and  then  interchanging  the  orders 
of  summation  and  expectation  in  the  last  two  terms,  we  get 

MSE  =  E  {d\n)]    -  2  Yk^)  E  {d{n)u{n  -  ^  +  1)} 

MM  ^^  (2.6) 

+    y  ^A(/t)A(m)  E  {u{n  -k+  l)u{n  -  m  -f-  1)} . 

Assuming  that  the  input  signal  u(n)  and  the  desired  response  d(n)  are  jointly  stationary, 
the  three  terms  on  the  right-hand  side  of  the  above  equation  may  be  interpreted  as  fol- 
lows: 

1.  The  expectation  E  {d^{n)}  is  equal  to  the  mean  square  value  of  the  desired  response 
d(n): 

P,=  E{d\n)].  (2.7) 

2.  The  expectation  E  {d{n)u{n  —  k  +  1)}  is  equal  to  the  cross-correlation  function  of 
the  desired  response  d(n)  and  the  input  signal  u(n)  for  the  lag  of  k-1.  We  can 
therefore  write  the  single  summation  term  on  the  righthand  side  of  Equation  (2.6) 
as  follows: 

M  M 

yh{k)E{d{n)u{n  -  ^  +  1)}  =  yh{k)p{k  -  1) .  (2.8) 

3.  Finally,  the  expectation  E  {u{n  —  k -\- \)u{n  —  m  +  I)]  is  equal  to  the 
autocorrelation  function  of  the  input  signal  u(n)  for  the  lag  of  m-k: 

r{m  -k)  =  E  {u{n  -k-\-  \)u{n  -m+\)}.  (2.9) 

Accordingly,  we  can  rewrite  the  double  summation  term  on  the  righthand  side  of 
Equation  (2.6)  in  the  form 

M        M 

y    Yjh{k)h{m)  E{u{n  -k+  l)u{n  -m+\)] 
^    ^"=1         MM  (2.10) 

=  y   ^/i(;t)/i(m)r(m-^). 


Thus,  substituting  Equations  (2.7),  (2.8),  and  (2.10)  into  Equation  (2.6),  we  find  that  the 
expression  for  the  mean  square  error  may  by  rewritten  in  the  form 

MSE=    P^   -  2  yh{k)p{k  -  I)  +    Y^  Yj^{k)h{m)r{ni  -  k) .  (2.11) 


k=i 


k^\    m=l 


By  differentiating  Equation  (2.11)  with  respect  to  h(k)  and  setting  it  equal  to  zero,  we 
have  the  following  set  of  M  simultaneous  equations 


p{k  -  1)  =  ^/io(m)r(m  -  k),      k=  1,2,...,M. 


(2.12) 


m=\ 


This  system  of  M  simultaneous  equations  is  called  the  normal  equations  with  optimum 
filter  coefficients  as  the  unknowns.   With  the  following  definitions, 
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(2.15) 


we  can  rewrite  the  normal  equations  of  Equation  (2.12)  in  matrix  form  as 


P  =  Rho- 


(2.16) 


This  represents  the  discrete-time  version  of  the  Wiener-Hopf  equation. 

B.     TWO  DIMENSIONAL  WIENER  FILTER 

The  following  derivation  parallels  the  development  by  Hadhoud  and  Thomas,  how- 
ever this  research  and  simulation  were  completed  separately  and  prior  to  the  publication 
of  Reference  9.  In  order  to  be  apphcable  for  an  image  processing  problem,  we  must 
extend  the  formulation  in  the  previous  section  to  two  dimensions  [Refs.  10,11].  This  is 
accomplished  by  developing  a  basic  two  dimensional  Wiener  filter  as  sho'wn  in 
Figure  2  on  page  9.  Within  this  filter,  we  use  two  input  images,  the  reference  array  X 
and  the  primary  input  array  D.  The  primary  input  image  D  is  a  two  dimensional  array 
which  represents  the  ideal  image  plus  additive  noise,  while  the  reference  image  X  is  noise 
that  is  assumed  to  be  correlated  to  the  noise  in  the  primary  input.  Both  the  input  arrays 
are  MxM  in  dimension.  The  Wiener  filter  is  an  NxN  causal  FIR  filter  with  a  set  of 
weights  Wj  defined  as 


Wj  = 


Wjil,0) 


Wj{0,\) 
^;(1,1) 


^/A^-1,0)    Wj{N- 1,1) 


Wj{0,N-\) 
Wj{l,n-l) 


Wj{N-l,N-l) 


(2.17) 


which  minimize  the  mean  of  the  squared  error,  e,,  between  the  filter  output  and  the  de- 
sired input  D.  We  designate  y  as  the  iteration  number  given  by  J  =  mM  +  n  where  m  and 
n  take  on  the  values  from  0  to  M-1.  Just  as  in  the  one  dimensional  case,  the  filter  out- 
put, y{m,n),  is  the  convolution  sum  of  the  filter  mask  and  the  reference  input  X  which  is 
given  by 


/V-l    N-\ 


yim,n)  =J]   Y  Wj{l,k)X{m  -  l,n  -  k). 


(2.18) 


/=0     k^ 


During  ihejth  iteration  the  input  from  array  X  is  represented  by  Xj  where 
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Figure  2.      T>vo  Dimensional  Wiener  Filter 


A}  = 


X{m,n)  X{m.n  —  1) 


X{ni  -  iV  +  1  ,n)    X{m  -  N+\,n-\) 


Xim,n-N+  1) 


X{m-N+l,n-N+  1) 


(2.19) 


Using  the  Equations  (2.17)  and  (2.19),  Equation  (2.18)  can  now  be  written  as 


/V-1    /V-1 


y{m,n)=y   YjWj{l,k)  Xj{l,k)  (2.20) 

for  ihtjth  iteration. 

This  output  y{m,n)  can  now  be  subtracted  from  one  pixel  D(m,n)  of  the  array  D  to 
produce  the  error  signal  at  ihe  Jih  iteration 

N-\    yV-1 

ej  =  D{m,n)  -  X   X  Wj{l,k)X{m  -  l,n  -  k).  (2.21) 

Since  the  purpose  of  the  Wiener  filter  is  to  provide  a  set  of  weights  which  minimize 
the  MSE,  we  can  denote  it  as 

MSE=^E{ej]  (2.22) 

where  £{.}  is  the  expectation  operator. 

Using  a  mathematical  derivation  similar  to  the  one  dimensional  case  of  the  previous 
section,  we  can  see  that 


N-\    N-\ 

ej  =  Z)^(w,«)  -  2  y    X  Wj{l,k)  D{m,n)  X{m  -  l,n  -  k) 

^r.l    N-\  N-\  N-i        i^  ^^=0  (2.23) 

+  Z  Z  Z  Z  ^^/^'^)  ^>'^^  ^(^  -  ^'"  -  ^^  ^(^  -  ^'"  -  ^^  • 

/=0       k=0    p=0     9=0 

Substituting  Equation  (2.23)  into  Equation  (2.22)  yields 


A-l    A^-1 

MSE  =  ElD\m,n)-]  -  2  Y    V  ^y(/,A)  £[Z)(m,rt)  X{m  -  l,n  -  A)] 

N-\    N-\    N-\    N-\  ^     ^  (2.24) 

+  Z  Z  Z  Z  ^/^'^)  ^^>'^^  ^^(^  ~  ^'"  ~  ^^  ^(^  ~  P'"  ~  ^^^  • 

Defining  P  as  the  crosscorrelation  matrix  between  the  desired  response  D(m,n)  and  the 
reference  input,  R  as  the  input  autocorrelation  matrix,  and  W  as  the  optimum  Wiener 
weight  matrix,  we  can  rewrite  Equation  (2.24)  as 
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/V-1    ,V-1 


^-l   ^'-l  A'-i  /V-i 


+ 


^    ^  ff^/a)  ^f^>,^)  i^(;.  -  l,q  -  k)  . 


7^     S^    ;>=0     9=0 


(2.25) 


Finally,  if  we  minimize  the  MSE  with  respect  to  Wj{l,k)  then  we  have 


A^-l    jV-1 


^(/''^)=Z  Z^>'^^^(^-^'^-^)- 


p=0     9=0 


(2.26) 


This  equation  is  the  two  dimensional  equivalent  to  Equation  (2.12).    The  matrix  form 
of  Equation  (2.26)  is  given  as 


P  =  RW 


(2.27) 


where  the  elements  of  this  equation  are  defined  as  follows 


R  = 


[^ol 


iRi'] 
\.R{] 


[^-yV+l]     [^-A'+2] 


[/?o] 


(2.28) 


Wo 


w  = 


(2.29) 


W 


N-] 


Px 


p  = 


(2.30) 


Pn-\ 
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Within  the  R  matrix  in  equation  (2.28)  each  element  is  a  block  Toeplitz  matrix  repres- 
ented by  the  equation 

R[p  -  l,q  -k)  =  ElX{m  -  l,n  -  k)X{m-p,n  -  ^)]  (2.31) 

where  l,k,p,q  range  from  0  to  N-1. 

C.     TWO  DIMENSIONAL  LEAST  MEAN  SQUARE 

One  means  of  obtaining  an  approximate  solution  for  the  optimum  weight  matrix, 
W,  is  the  use  of  a  two  dimensional  LMS  algorithm  which  is  depicted  in  Figure  3  on  page 
13.  This  differs  from  Figure  2  on  page  9  in  that  the  error  ej  is  used  to  update  the  filter 
coefficients  before  shifting  the  data  window  Xj  across  the  reference  input  for  the  next 
iteration.  As  in  the  one  dimensional  LMS  algorithm,  we  are  updating  the  coefficient 
matrix  by  adding  the  present  matrix  to  a  change  proportional  to  the  negative  gradient 
of  the  error  where  the  one  dimensional  instantaneous  estimates  of  the  gradient  vector 
are  based  on  sample  values  of  the  input  and  the  error  signal  e{n)u{n)  [Ref  8].  For  the 
two  dimensional  Jf/i  iteration,  we  define  the  updated  matrix  as 

Wj^,  =r  iVj  -  fx  ej  Xj  (2.32) 

where 

Wj^i  updated  weight  matrix 

Wj  previous  weight  matrix 

IX  scaler  multiplier  controlling  the  rate  of  convergence  and  filter  stability 

{ej){Xj)  estimate  for  the  2-D  instantaneous  gradient 

The  previous  equation  can  also  be  written  as 

Wj^^{l,k)  =  Wj{l,k)  +  2  ^  (ej)  X{m  -  l,n  -  k).  (2.33) 

These  two  equations  give  the  two  dimensional  weight  updating  algorithm  for  the  LMS 
adaptive  filter.  The  algorithm  we  have  developed  here  may  be  implemented  "without  any 
form  of  matrix  operations,  averaging,  or  differentiation. 

The  value  of  /z  may  be  chosen  based  on  the  desired  tracking  ability,  steady-state 
mean  square  error,  and  convergence  speed.  In  signal  processing,  there  are  several 
methods  for  determining  a  suitable  value,  however  in  many  of  these  cases  it  requires 
knowledge  of  the  eigenvalues  and  eigenvectors  of  the  autocorrleation  matrix.    In  image 
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Figure  3.      T»o  Dimensional  Least  Mean  Square 

processing,  one  method  which  does  not  require  these  values  for  the  computation  of/j  is 
trial  and  error  based  on  output  image. 

An  alternate  method  used  in  one  dimensional  design  which  again  does  not  require 
a  priori  knowledge  of  the  autocorrelation  matrix,  discussed  by  Bitmead  and  Anderson 
[Ref  12],  is  the  normalized  LMS.  Using  our  previous  notation,  this  method  can  be  ex- 
tended to  two  dimensions  by  first  considering  the  two  dimensional  LMS  update  equation 
(2.32).    In  this  equation,  we  redefine  the  step-size  parameter,  ix,  for  a  given  /  and  k  as 
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m(/./c)  = 


a + o]{i,k) 


(2.34) 


in  which  a,  represents  the  normalized  step  size  chosen  between  zero  and  two,  /?  is  an- 
other small  positive  constant,  and  a]{l,k)  is  one  of  the  values  from  the  matrix 


aj  = 


a]{0,0)  a]{0,\) 

a](l,0)  ^^(l.l) 


o](0,N-l) 
oj{l,N-l) 


a]{N-\,0)    ajiN-l,l)    ...    a]{N -  l,N -  I) 


(2.35) 


The  element  required  from  the  matrix  depends  upon  the  current  values  of  /  and  k  being 
used  for  equation  (2.33).  The  matrix  may  be  be  intitialized  with  the  values  of 
{X{m  —  l,n  —  k)Y  and  to  update  the  values  in  a],  we  use  the  following  equation 


<7]+i(a)  =  P{o]{l,k))  +  {\-p){X{m-l,n-k)f 


(2.36) 


where  p  is  a  weighting  factor  between  zero  and  one.    This  ensures  that  the  value  of  m 
does  not  become  large  enough  to  cause  the  algorithm  to  become  unstable. 

D.     IMPLEMENTATION 
1.     System  Identification 

In  order  to  initially  test  the  two  dimensional  least  mean  square  algorithm,  we 
established  a  system  identification  model  shown  in  Figure  4  on  page  15.  Within  this 
model  the  output  of  the  known  FIR  filter,  d{m,n),  is  defined  as 


d{m,n)  =  .4  W{m,n)  +  .6  W{m-  l,n) 
-  .3  W{m,n  -  1)  +  .2  W{m  -  l,n  -  1) . 


(2.37) 


The  output  of  the  two  dimensional  least  mean  square  filter,  y{m,n)  ,  consists  of  a  set  of 
adjustable  coefficients  and  is  defined  as 


y{m,n)  =  AO  W{m,n)  +  Al  W{m-  \,n) 
+  A2  W{m,n  -  \)  +  Al  W{m  -  l,n  -  1) . 


(2.38) 


The  adaptive  filter  output  y{m,n)  is  then  compared  with  the  known  system  output 
d{m,n)  to  produce  an  error  signal  e{m,n),  defined  as  the  diiference  between  them. 
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Figure  4.      T>\o  Dimensional  System  Identification 


The  operation  of  the  adaptive  filter  is  to  minimize  the  error  signal  e{m,n)  by 
providing  an  adaptive  process  for  adjusting  the  coefficients,  for  this  adaptive  process 
we  use  the  update  equation  (2.32)  for  the  two  dimensional  least  mean  square  algorithm 
developed  in  the  previous  section 


W, 


J+i 


IV 


nejXj 


(2.32) 


For  this  model  a  32x32  white  gaussian  noise  matrix  was  used  as  the  input  and 
the  rate  of  convergence  can  be  be  seen  in  Figure  5  on  page  16  and  Figure  6  on  page 
17.  Following  600  iterations  all  the  coefficients  had  converged  to  within  10"^  of  the  ac- 
tual coefficient  value  and  the  error  was  3j[:10-^  A  computer  program  for  this  system 
identification  model  is  given  in  Appendix  A. 
2.     Adaptive  Noise  Canceler 

The  usual  method  of  estimating  a  signal  corrupted  by  additive  noise  is  to  pass 
the  composite  signal  through  a  filter  that  tends  to  suppress  the  noise  while  leaving  the 
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Figure  5.      LMS  System  Identification-Rate  of  Convergence 

signal  relatively  unchanged.  The  noise  canceler  and  adaptive  hne  enhancer  developed 
by  Widrow  [Rcf.  13]  are  well-documented  ways  of  doing  that.  Adaptive  noise  canceling 
is  a  variation  of  optimal  filtering  that  is  highly  advantageous  in  many  applications.  It 
uses  an  auxiliary  or  reference  input  derived  from  one  or  more  sensors  located  at  points 
in  the  noise  field  where  the  signal  is  weak  or  undetectable.  This  input  is  filtered  and 
subtracted  from  a  primary  input  containing  both  signal  and  noise.  The  reference  input 
and  the  noise  in  tiie  primary  input  are  therefore  correlated,  and  as  a  result  the  primary 
noise  is  attenuated  or  eliminated  bv  cancellation. 


16 


?^ 


o 


o 


8, 

■  a 

o 
I 


o 
I 


LEC^ND 
ERROR 


— T— 

100 


I  I 

200       900 


400       600       600 

ITERATIONS 


700       000       900      1000 


Figure  6.      LMS  System  Identification-Rate  of  Convergence 

The  basic  noise  canceling  concept  is  illustrated  in  Figure  7  on  page  18.  A  signal 
is  transmitted  over  a  channel  to  a  sensor  that  receives  the  signal  plus  noise,  Hq.  The 
combined  signal  and  noise  s  +  n^  forms  the  "primary  input"  to  the  canceler.  A  second 
sensor  receives  a  noise  /j,  which  is  correlated  in  some  way  with  the  noise  /Iq.  This  sensor 
output  provides  the  "reference  input"  to  the  canceler.  The  noise  n^  is  filtered  to  produce 
an  output }'  that  is  a  close  replica  of  Wq  .  This  output  is  subtracted  from  the  primary 
input  5  +  Wq  to  produce  the  system  output  5  +  Hq  —y. 
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Figure  7.      Adaptive  Noise  Canceler 

In  the  system  shown  in  Figure  7,  the  reference  input  is  processed  by  our  two 
dimensional  least  mean  square  filter.  Thus  the  filter  operates  under  changing  conditions 
and  will  readjust  itself  continuously  to  minimize  the  error  signal. 

The  computer  program  which  simulates  this  noise  canceling  model  is  provided 
in  Appendix  B.  We  will  discuss  the  results  and  conclusions  concerning  various  simu- 
lations in  Chapter  4. 

3.     Adaptive  Line  Enhancer 

A  special  case  of  adaptive  noise  canceling  is  when  there  is  only  one  signal  x{n) 
available  which  is  contaminated  by  noise.  In  such  a  case,  the  signal  x[j\)  provides  its 
own  reference  signal  >(/?),  which  is  taken  to  be  a  delayed  replica  of 
x{n):y{n)  =  x{n  —  A),  as  shown  in  Figure  8  on  page  19.  The  adaptive  filter  will  respond 
by  canceling  any  components  of  the  main  signal  x{n)  that  are  in  any  way  correlated  with 
the  secondary  signal  y{n)  =  x{n  —  A).  Suppose  the  signal  x{n)  consists  of  two  compo- 
nents: a  narrowband  component  that  has  long-range  correlations  and  a  broadband 
component  which  will  tend  to  have  short-range  correlations.  One  of  these  could  repre- 
sent the  desired  signal  and  the  other  an  undesired  interfering  noise.    Suppose  that  the 
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Figure  8.      Adaptive  Line  Enhancer 

delay  A  is  selected  so  that  it  falls  between  the  correlation  lengths.  Since  A  is  longer  than 
the  elVective  correlation  length  of  the  broadband  component,  the  delayed  replica  will  be 
entirely  uncorrected  with  the  broadband  part  of  the  main  signal.  The  adaptive  filter 
will  not  be  able  to  respond  to  this  component.  On  the  other  hand,  since  A  is  shorter 
than  the  correlation  length  of  the  narrowband  component,  the  delayed  replica  that  ap- 
pears in  the  secondary  input  will  be  correlated  with  narrowband  part  of  the  main  signal 
and  the  filter  will  respond  to  cancel  it.  Note  that  if  A  is  selected  to  be  longer  than  both 
correlation  lengths,  the  secondary  input  will  become  uncorrected  with  the  primary  input 
and  the  adaptive  filter  will  turn  itself  off.  In  the  opposite  case,  when  the  delay  A  is  se- 
lected to  be  less  than  both  correlation  lengths,  then  both  components  of  the  secondary 
signal  will  be  correlated  with  the  primary  signal,  and  therefore  the  adaptive  filter  will 
respond  to  cancel  the  primar>'  x{n)  completely.  The  computational  algorithm  for  the 
adaptive  line  enhancer  is  shown  in  the  following  three  equations: 


\f 


M 


x{n)  =  2^/7^(«)v(/7  -  m)  =  2_^hjn).x{n 


m 


A) 


(2.40) 


m=0 


wi=0 
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e(n)  =  x{n)  -  x{fi)  (2.41) 

hJn+l)  =  hJn)  +  2^iein)x{n-m-A)       m  =  0,1,2,.. .,M  (2.42) 

For  this  model  we  also  developed  a  computer  simulation  which  incorporates  our 
two  dimensional  LMS  algorithm  and  it  is  provided  in  Appendix  C.  The  results  and 
conclusions  will  again  be  discussed  in  Chapter  4. 
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III.      TWO  DIMENSIONAL  ADAPTIVE  RECURSIVE  LEAST  SQUARES 

ALGORITHM 

In  the  previous  chapter,  under  the  LMS  algorithm,  the  available  data  samples  were 
used  in  order  to  attempt  to  move  the  current  estimate  of  the  impulse  response  to  the 
optimum  value.  This  approach  has  the  advantage  of  being  simple  to  implement  but 
carries  with  it  the  disadvantages  that  it  can  be  slow  to  approach  the  optimal  weight 
vector  and,  once  close  to  it,  will  usually  fluctuate  around  the  optimal  vector  rather  than 
actually  converge  to  it  due  to  the  effects  of  approximations  made  in  the  estimate  of  the 
performance  function  gradient. 

To  overcome  these  difficulties,  we  examine  another  approach  in  this  chapter  which 
uses  the  input  data  in  such  a  way  as  to  ensure  optimality  at  each  step.  This  alternative 
algorithm  is  based  on  the  exact  minimization  of  the  least  square  criteria.  This  algorithm 
is  known  as  recursive  least  squares  (RLS). 

A.     ONE  DIMENSIONAL  RECURSIVE  LEAST  SQUARES 

As  in  the  LMS  case,  the  one  dimensional  RLS  algorithms  is  developed  in  several 
different  ways  [Refs.  5,  8  ].  Orfanidis  [Ref  6]  provides  a  derivation  which  we  will  con- 
sider prior  to  our  extension  to  two  dimensions.  The  tapped  delay  line  shown  in 
Figure  9  on  page  22  will  provide  the  reference  for  the  following  discussion.  We  begin 
by  replacing  the  LMS  estimation  criteria  of  MSE  =  £[e^]  by 


n 


MSE  =    )  e\k)        A  =  0,  l,...,n  (3.2) 


k=0 

where 


e{k)  =  x(k)  -  xik)  (3.3) 

and  x{k)  is  the  estimate  of  jr(/c)  produced  by  the  Mth-order  Wiener  filter 


.M 


X  =  2^h{m)y{k-m)  (3.4) 

m=0 
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Figure  9.      One  Dimensional  RLS  Reference  (TDL) 

Substituting  equation  (3.2)  into  (3.1)  and  setting  the  derivative  with  respect  to  h  to  zero 
we  find  the  least-squares  analogs  of  the  orthogonality  equations 
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n 

dMSE 2ye{k)y{k)  =  0  (3.5) 


^h  ;^o 


which  may  be  rewritten  in  their  normal  equation  form  as  follows 


n 

I 


{x{k)  -  h''y{kMk)  =  0  (3.6) 


r  "  "I  " 

^y(/:)y(/c)^    h  =  yx{k)y{k)  (3.7) 


Defining  the  quantities 


R{")  =  }  yii<)yikf  (3.8) 


^^ 


n 


r(«)  =    }x{k)  y{k)  (3.9) 

we  can  then  write  the  normal  equation  as  R{n)  h  =  r{n),  with  solution  h  =  R{n)-^  r{n). 
Note  that  the  n-dependence  of  R{n)  and  r{n)  depend  on  the  current  time  n,  therefore, 

h(/7)  =  R{n)~\{n)  =  P{n)r{n)  (3.10) 

where  P{n)  =  R{n)-\  These  are  the  least  squares  analogs  of  the  ordinary  Wiener  sol- 
ution, with  R{n)  and  r(A7)  playing  the  role  of  the  covariance  matrix  •£Iy(«)y(«)T  and 
cross-correlation  vector  Elx{n)y{n)'],  respectively.  The  RLS  algorithm  is  obtained  by 
writing  equation  (3.10)  recursively  in  n  and  then  using  the  following  matrix  inversion 
lemma  [Ref.  5] 

{A  +  BCD)'^  =  A~^  -  A~^B{DA~^B  +  C'^Y^DA'^  (3.11) 

we  get  the  update  equation  for  the  P  matrix 

Pin)  ^  Pin -X)-     ^(-Oy(;)y(-)^/-(---l).  (3,2) 

1  +  y{n)'P{n-\)y{n) 
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Using  the  quantities  in  equations  (3.8)  and  (3.9)  to  satisfy  the  recursions  yields 

R{n)  =.  R{n-\)  +  y{n)y{nf  (3.13) 

r{n)  =  r(«-l)  +  x{n)y{n)  (3.14) 

and  substituting  equations  (3.12)  and  (3.14)  into  (3.10)  after  some  mathematical  manip- 
ulations we  fmd 

h{n)  =  h(«  -  1)  +  P{n)  y{n)  e{n)  (3.15) 

which  differs  from  the  LMS  alogorithm  by  the  presence  of  the  P(n),  vice  /x,  in  front  of 
the  weight  correction  term.  Since  P{n)  =  R{n)-^  is  a  measure  of  the  covariance  matrix 
E[y{n)y{n)^,  the  presence  of  R{n)~^  makes  the  RLS  algorithm  behave  like  Newton's 
method,  and  hence  has  fast  convergence  properties. 

B.     TWO  DIMENSIONAL  RECURSIVE  LEAST  SQUARES 

Although  it  is  discussed  briefly  by  Wellstead  and  Caldas  Pintos  [Ref  14],  limited 
information  is  available  in  the  open  literature  in  this  area.  In  order  to  develop  a  two 
dimensional  recursive  least  square  (RLS)  algorithm  we  will  extend  the  one  dimensional 
adaptive  algorithm  developed  in  the  previous  section  in  a  method  similar  to  that  used  for 
the  LMS  algorithm. 

As  in  Chapter  2,  using  Figure  10  on  page  25  as  a  reference  we  see  that  the  basic 
filter  has  two  input  images.  The  primary  two  dimensional  input  array,  D,  is  the  ideal 
image  plus  additive  noise.  The  reference  image  X  is  the  noise  array.  Each  array  is  of 
dimension  M  by  M.   The  filter  mask,  W,  is  N  by  N. 

The  one  dimensional  RLS  algorithm,  equation  (3.15),  is  the  same  as  the  one  di- 
mensional LMS  algorithm  with  the  exception  that  P{n)  replaces  n.  Therefore,  if  a 
method  of  developing  a  P  matrix  for  the  two  dimensional  case  can  be  devised  we  can 
then  use  an  algorithm  similar  to  the  two  dimensional  LMS  algorithm  to  update  the  filter 
coefficients.   First,  let 
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Figure  10.      T>vo  Dimensional  RLS  Reference 
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Pj- 


P/0,0)  Pj{0,\) 


Pj{^,N^-\) 


Pj{N^ -1,0)   PjiN'-l,\)    -pj{N^-l,N^-l) 


(3.16) 


which  represents  the  P  matrix  on  the  yth  iteration  where  J  is  defined  as  j  =  mM  +  n. 
As  in  the  two  dimensional  LMS  algorithm,  we  let  the  filter  coefficient  matrix  be  given 
as 


Wj  = 


Wj{0,0) 
Wj{\,0) 


Wj{0,\) 

Wj{\,\) 


Wj{N-\,0)    Wj{N-\,\) 


Wp,N-\) 
Wj{\,N-\) 


Wj{N-\,N-\) 


(3.17) 


and  the  input  data  window  as 


A}  = 


X{m,n) 


X{m,n  -  1) 


X{m-N^\,n)    X{m-N+\,n-\) 


X{m,n-N-^\) 


X{m-N+\,n-N+\) 


(3.18) 


In  order  to  establish  an  update  equation  for  the  P  matrix  in  which  each  element  will 
have  the  proper  dimensions,  we  use  equations  (3.17)  and  (3.18)  and  make  the  following 
transformations;  the  filter  mask,  Wj,  equation  (3.17)  is  transformed  into  a  vector  defined 
as 


26 


IVWj  = 


Wj{\,N-\) 


lVj{N-l,\) 


Wj{N-l,N-l) 


and  we  transform  the  Xj  matrix,  equation  (3.18)  into 


X{m,n) 
X{m,n  -  1) 

X{m,n-N+  1) 


XXj  = 


X{m-N+\,n) 
X{m-N-^  1,«-  1) 

X{m-N+  \,n-N+\) 


Using  the  quanitities  defmed  above,  we  can  then  write  the  P  matrix  update  equation  as 


(3.19) 


(3.20) 
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Pj  =  PJ-^  - 


Pj_,{xxj){xx;)Pj_, 

1  +  {XXj)Pj_,{XXj) 


(3.21) 


As  discussed  earlier,  the  one  dimensional  RLS  weight  updating  algorithm  difiers 
from  the  LMS  in  that  it  contains  the  P  matrix  while  the  LMS  algorithm  contains  n.  This 
must  be  considered  in  two  dimensions  since  /*  is  a  matrix  and  ^  is  a  scalar.  Therefore 
using  the  previous  equations,  we  can  obtain  the  update  to  the  filter  coefficients  via 


WWj^,  =   WWj  -  {Pj){ej){XXj) 


(3.22) 


where  ej  is  given  as 


ej  =  D{m,n)  -  y{m,n) 


(3.23) 


y{m,n)  is  the  convolution  sum  of  the  image  in  the  reference  input  with  the  filter  window 


A-l      A^-1 


(3.24) 


and  the  value  P,  is  updated  by  equation  (3.21). 

The  cost  function  of  the  RLS  criterion  could  be  modified  to  include  a  windowing 
of  the  input  data  as  follows 


MSE  = 


n 

k=0 


e  (k) 


This  modification  results  in  the  following  modified  P  matrix  update: 


P  --L 


^y-i- 


Pj_,{xxj){xx;)Pj_, 

\ -^  {XXj)Pj_,{XXj) 


(3.25) 


(3.26) 


where  q,  the  averaging  or  "forgetting  factor"  is  a  positive  constant.  It  is  usually  chosen 
to  be  slightly  less  than  one,  thereby  diminishing  the  contribution  of  the  "older"  data. 
The  problem  of  data  nonstationarity  is  the  main  reason  for  introducing  this  type  of 
weighting  factor.[Ref  6]  It  should  be  noted  that  the  usual  RLS  algorithm  is  attained 
when  q  =  1  and  that  no  changes  in  the  amount  of  computation  is  required. 
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C.     IMPLEMENTATION 

1.     System  Identification 

Following  the  mathematical  development  of  the  two  dimensional  recursive  least 
squares  algorithm,  we  implemented  and  tested  it  as  was  previously  done  with  the  two 
dimensional  least  mean  square  algorithm.  Our  initial  testing  was  done  in  a  system 
identification  model  as  shown  in  Figure  11.  As  in  the  two  dimensional  LMS  case,  we 
used  a  known  filter  with  the  following  output  equation 


^(,„,„)  =  .4lV{m,n)  +   .6JF(m-  1,«) 
-  3lV{m,n-\)  +  .2lV{m-  \,n-  1) 


(3.27) 


WHITE 

GAUSSIAN 

NOISE 


KNOWN 
FILTER 

d(ni,n) 

^ 

ADAPTIVE 
RLS  FILTER 


~Z 


y(ni,n) 


error 


Figure  11. 


Two  Dimensional  RLS  System  Identification 
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The  two  dimensional  adaptive  recursive  least  square  filter  output  was  given  as 

y{m,n)  =  AO  W{m,n)  +  A I  W{m-l,n) 

+  A2  W{m,n  -  1)  +  ^3  W{m  -\,n-\)  ^^'^^^ 

The  error  signal  was  generated  by  taking  the  difference  between  the  desired  (known) 
signal,  d{m,n),  and  the  filter  output,  j;(m,rt).  By  using  equations  (3.21)  and  (3.22)  to  up- 
date the  P  matrix  and  the  filter  coefficients,  we  caused  the  filter  weights  to  move  toward 
their  optimum  values  and  the  error  signal  to  approach  zero. 

The  computer  program  for  this  model  is  provided  is  Appendix  D.  For  this  sys- 
tem, a  32x32  white  gaussian  noise  matrix  was  generated  for  the  input  to  both  the  known 
filter  and  the  RLS  adaptive  filter.  The  rate  of  convergence  is  shown  in  Figure  12  on 
page  31  and  Figure  13  on  page  32  and  as  can  be  seen  after  approximately  70  iterations 
each  of  the  coefficients  had  converged  to  with  10"^  of  the  actual  value  and  the  error  was 
10-^ 

2.     Adaptive  Noise  Canceler/ Adaptive  Line  Enhancer 

These  two  systems  were  both  discussed  in  Chapter  2  with  regard  to  implemen- 
tation of  the  two  dimensional  LMS  algorithm.  In  order  to  further  test  the  RLS  algo- 
rithm we  also  implemented  it  within  the  noise  canceler  and  the  adaptive  line  enhancer. 
The  computer  programs  which  were  used  for  the  simulation  are  given  Appendix  E  and 
Appendix  F,  respectively.  The  simulation  results  and  discussion  will  be  provided  in 
Chapter  4. 
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Figure  12.      RLS  System  Identification  Rate  of  Convergence 
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IV.     RESULTS  AND  CONCLUSIONS 

In  this  chapter,  we  present  the  experimental  results  from  our  implenientation  of  the 
two  dimeiisionai  LMS  algorithm  derived  in  Chapter  II  and  the  two  dimensional  RLS 
algorithm  derived  in  Chapter  HI.  The  two  algorithms  were  used  in  a  noise  canceler  and 
an  adaptive  line  enliancer,  as  discussed  in  Chapter  II,  to  solve  an  image  restoration 
problem. 

Figure  14  shows  the  original  image  of  a  house  which  is  comprised  of  128  by  128 

eight  bit  pixels.    The  image  has  a  mean  value  of  131.68  and  a  variance  of  3194.4.    Li 

Figure  15  on  page  34.  we  have  corrupted  the  original  image  with  additive  white  gaiissian. 

noise  (zero  mean  and  variance  1600). 

"1 


igure   14,      Original  Image 
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rf9'r''tutfli^ffiv^rruimi!pr* 


Figure  15.      image  plus  Additive  Noise 

.4.     NOISE  CANCELER  RESULTS 

In  order  to  soive  our  image  restoration  problem,  we  initially  implemented  botli  the 
LMS  and  the  RLS  algoritiims  in  a  noise  canceler  using  a  two  by  two  filter  matrix  and 
a  value  of  70a- iO"-  for  ij.  in  the  LMS  algorithm, 

Figure  16  on  page  35  shows  the  output  when  the  LMS  algorithm  is  processed 
through  the  image  for  one  pass;  this  results  in  16,384  updates  of  the  filter  corresponding 
10  the  16,384  pixels.  Figure  17  on  page  36  is  the  results  following  two  passes  through 
the  image.   No  significant  improvement  was  realized  with  more  than  two  passes, 

Lmpletneniing  a  two  by  two  R.LS  adaptive  filter  in  the  noise  canceler  produced  very 
favorable  results  after  one  pass  through  the  imiage.  These  results  are  shown  in 
Figure  IS  oti  page  36  and  compare  well  with  the  original  image  and  the  LMS  output 
after  two  passes. 

For  the  LMS  algorithm  our  best  results  were  achieved  with  a  value  of  lOxlQ-^  for 
fx.  however  in  order  to  demonstrate  the  effect  of  changing  this  value  Figure  19  on  page 
37  represents  the  output  v,'iih  ,a  ~  35.y10-'  and  Figure  20  on  page  37  is  produced  by  a 
spatially  varying  normalized  f.i  as  discussed  in  Chapter  IF   Ii  can  be  seen  that  for  various 
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values  of  a,  dilTerent  features  within  tlie  image  were  restored  at  diilerent  levels.  In  Fig- 
ure 16  on  page  35,  the  edges  and  more  detailed  segments  of  the  image  were  restored 
better  than  in  Figure  19  on  page  37,  however  areas  of  similar  contrast  v/ere  not  as  well 
restored.  When  using  tho,  normalized  fi  the  mean  value  more  closely  approached  the 
original  mean,  however  the  variance  v/as  reduced. 

Finally,  increasing  the  number  of  coefficients,  using  a  three  by  three  filter  matrix  in 
the  LMS  algorithm  we  obtained  the  results  shown  in  Figure  21  on  page  38  in  one  pass. 
This  filter  showed  no  significant  improvement  in  the  variance  compared  to  the  two  by 
two  filter;  however,  it  did  improve  with  respect  to  the  mean. 

Table  1  on  page  38  shows  the  restoration  results  comparing  the  mean  and  variance 
of  the  various  outputs  with  that  of  the  original  image  and  the  image  plus  noise. 


Figure  16,      Restored  Image:  Noise  Canceier/ 


^ne  pass) 


35 


Figure   17.      Restored  Image:  Noise  Cauceler/LMS  (Two  passes) 


Figure   IS,      Restored  Image:  Noise  Caiice!er/RLS  (One  pass) 
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Figure  19=      Restored  image:  Noise  Canceier/LMS  (different  mu) 


Figure  20.      Restored  Image:  Noise  Canceler/LMS  (normalized  mii) 
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Figure  21.      Restored  Image:  Noise  Canceler/LMS  (3  by  3) 


Table  1.     NOISE  CANCELER  RESULTS 


TYPE  OF  FILTER 

MEAN 

VARIANCE 

Original  Image 

131.68 

3194.4 

Image  plus  Noise 

134.25 

1647.7 

One  pass  2x2  LMS 

114.52 

2034.3 

Two  pass  2x2  LMS 

136.78 

3711.8 

One  pass  2x2  RLS 

139.04 

3607.7 

One  pass  L\lS(difrerent 
mu) 

120.12 

1940.4 

One  pass 

LMS( normalized  mu) 

126.53 

1500.7 

One  pass  3x3  LMS 

119.13 

2017.3 
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B.     ADAPXn'E  LINE  ENHANCER  RESULTS 

We  next  implenienied  our  two  algorithms  in  an  adaptive  line  enhancer.  In  order  to 
obtain  a  reference  input  from  the  primary  input,  we  used  a  two  dimensional  delay  oper- 
ator or(ni-i.n).  As  in  the  noise  canceler  we  again  used  a  tv/o  by  two  filter  matrix  and 
this  produced  the  results  shown  in  Figure  22  on  page  39  and  Figure  23  on  page  40  for 
the  one  pass  LMS  and  the  one  pass  RLS,  respectively.  No  significant  improvement  was 
noted  in  the  two  pass  LMS  compared  to  the  one  pass  LMS, 

Table  2  on  page  40  shows  the  comparison  of  each  output  mean  and  variance  with 
the  mean  and  variance  of  the  original  image  and  the  image  plus  noise. 


Figure  22.      Restored  Image:  Adaptive  Line  Enhancer/ LMS  (One  pass) 


L 


Figure  23.      Restored  Image:  Adaptive  Liiie  Enhancer/ RLS  (One  pass) 


Table  2.     ADAPTIVE  LINE  ENHANCER  RESULTS 
iPE  OF  FILTER 


Oriiiinal  linage 


linage  plus  Noise 


One  pass  2x2  LMS 


Two  pass  2x2  LMS 


I  One  pass  2x2  R.LS 


C.     CONCLUSIONS 

ill  this  thesis  we  have  extended  a  one  dimetisionai  LMS  and  a  one  dimensional  RLS 
adaptive  algorithm  to  tv^'o  dimensions.  As  we  examine  the  results  we  see  that  many  of 
the  one  dimeDsional  comparisons  between  LMS  and  RLS  are  applicable  to  two  dimen- 
sions; 

L  the  rate  of  convergence  of  the  RLS  a'gorithm  is  in  general  faster  than  that  of  the 
LMS  algorithm  by  an  order  of  magnitude 

2.  this  superior  performance  of  the  RLS  algorithm  is  attained  at  the  expense  of  a  large 
increase  in  the  computational  complexity 
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3.  there  are  no  approximations  within  the  derivation  of  the  RLS  algorithm,  therefore 
as  the  number  of  iterations  approach  infinity,  the  least  squares  estimate  approaches 
the  optimum  Wiener  value 

4.  in  the  RLS  algorithm  the  correction  which  is  applied  is  computed  using  all  past 
data  whereas  the  LMS  algorithm  uses  only  the  instantaneous  sample  and  the  error 
signal;  this  is  not  necessarily  an  advantage  in  image  processing. 

The  objectives  of  this  thesis  were  successfully  completed.  Some  suggestions  for  fu- 
ture work  include:  (i)  to  examine  these  two  dimensional  algorithms  in  other  applica- 
tions, applying  them  to  areas  other  than  image  processing,  (ii)  to  extend  the  two 
dimensional  LMS  and  RLS  algorithms  to  n-dimensions  and  implement  them,  (iii)  to 
analyze  the  extension  of  these  one  dimensional  algorithms  to  two  dimensions  in  the  fre- 
quency domain.  Multidimensional  digital  signal  processing  is  being  applied  in  many 
areas  today,  however  the  potential  for  future  growth  and  applications  appears  unlimited. 
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APPENDIX  A.     LMS  SYSTEM  IDENTIFICATION  PROGRAM 

C     SYSTEM  IDENTIFICATION  OF  AN  ADAPTIVE  TWO  DIMENSIONAL 
C         LMS  FILTER  VS.  A  KNOWN  TWO  DIMENSIONAL  FILTER 
C 

REAL  MU 

DIMENSION  A0(0: 32,0:  32) ,A1(0:  32,0:  32) ,A2(0:  32,0:  32) , 

*  A3(0:  32,0:  32)  ,E(0:  32,0:  32)  ,D(0:  32,0:  32)  ,W(-1:  32,-1:  32) , 

*  Y(0:32,0:32) 
C 

C     ******VARIABLE  DEFINITIONS****** 
C     A=NOISE  AMPLITUDE 
C     S=NOISE  VARIANCE 
C     AM=NOISE  MEAN 
C     IX=SEED 

C     A0,A1,A2,A3=FILTER  COEFFICIENTS 
C     MU=SCALING  FACTOR  WITHIN  THE  LMS  ALGORITHM 
C     W=WHITE  GAUSSIAN  NOISE  GENERATED  BY  SUBROUTINE 
C     Y=FILTER  OUTPUT 
C     D=DESIRED  OUTPUT 
C     E=ERROR(D  -  Y) 
C 

C     Vf*yfy-yf*yf INITIAL  VALUES******* 
A=l. 
MU=.  01 
A0(0,0)=0. 
A1(0,0)=0. 
A2(0,0)=0. 
A3(0,0)=0. 
S=l. 
AM=0. 
IX=65539 
C 

C     yc***ycOPEN  FILE  AND  FILL  INPUT  MATRIX  WITH  WHITE 
GAUSSIAN  NOISE******** 
OPEN  (UNIT=80,FILE='SYSID2  DATA ', STATUS= ' NEW ' ) 
DO  10  M=-l,31 
DO  20  N=-l,31 

IF((M.  LT.  0).0R.  (N.  LT.  0))THEN 

W(M,N)=0. 
ELSE 

CALL  GAUSS(IX,S,AM,V) 
W(M,N)=V*A 
ENDIF 
20       CONTINUE 
10    CONTINUE 
C 

C     ******COMPUTE  THE  UPDATED  FILTER  COEFFICIENTS 
C         USING  THE  TWO  DIMENSIONAL  LMS  ALGORITHM***** 
DO  40  K=0,31 
DO  50  J=0,31 

D(K,J)=. 4*W(K, J)+. 6*W(K-1, J)-. 3*W(K,J-1)+.  2*W(K-1,J-1) 
Y(K,J)=A0(K,J)*W(K,J)+A1(K,J)*W(K-1,J)+A2(K,J)*W(K,J-1)+ 
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*  A3(K,J)*W(K-1,J-1) 
E(K,J)=D(K,J)-Y(K,J) 
A0(K,J+1)=A0(K,J)+MU*E(K,J)*W(K,J) 
A1(K,J+1)=A1(K,J)+MU''^E(K,J)*W(K-1,J) 
A2(K,J+1)=A2(K,J)+MU'VE(K,J)''«W(K,J-1) 
A3(K,J+1)=A3(K,J)+MU*E(K,J)''^W(K-1,J-1) 

50      CONTINUE 

A0(K+1,0)=A0(K,J) 
A1(K+1,0)=A1(K,J) 
A2(K+1,0)=A2(K,J) 
A3(K+1,0)=A3(K,J) 
40  CONTINUE 
C 

Q     ******0UTPUT  RESULTS******* 
L=0 

DO  60  K=0,31 
DO  70  J=0,31 
L=L+1 

PRINT  15,  L,E(K,J),A0(K,J),A1(K,J),A2(K,J),A3(K,J) 
WRITE  (80,15)  L,E(K,J),A0(K,J),A1(K,J),A2(K,J),A3(K,J) 
15  FORMAT  ('  ' ,1X,I4,2X,F8.5,2X,F8.5,2X,F8.5,2X,F8.5, 

*  2X,F8.5) 
70       CONTINUE 

60    CONTINUE 

STOP 

END 
C 

SUBROUTINE  GAUSS(IX,S,AM,V) 

A=0.0 

DO  50  1=1,12 

CALL  RANDU(IX,IY,Y) 

IX=IY 
50  A=A+Y 

V=(A-6.0)*S+AM 

RETURN 

END 


C 
C 


SUBROUTINE  RANDU( IX, IY,YFL) 

IY=IX*65539 

IF(IY)5,6,6 

5  IY=IY+2147483647+1 

6  YFL=IY 

YFL=YFL*.  4656613E-9 

RETURN 

END 
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APPENDIX  B.     LMS  ALGORITHM  IMPLEMENTED  IN  A  NOISE 

CANCELER 

C     THIS  IS  A  VAX/ VMS  FORTRAN  PROGRAM  THAT  IMPLEMENTS  A  TWO 

C     DIMENSIONAL  2X2  ADAPTIVE  LMS  FILTER  WITHIN  A  NOISE 

C     CANCELER 

C 

INTEGER  M,N,K,J 

BYTE  B(0:  127) ,BYTEE(0:  127) ,BYTEU(0:  127), 

INTEGER*4  INTE(0:  127,0:  127),INTU(0:  127,0:  127), 

*  IE,IU 

REAL*4  MU,A,FMINE,FMAXE,FMINU,FMAXU 

REAL*4  AA(0:3),E(0:  127,0:  127), U(0:  127,0:  127) 

REAL*4  W(-l:  127,-1:  127), Y(0:  127,0:  127),IM(0:  127,0:  127) 
C 

C     ***.v^v*vARIABLE  DEFINITIONS****** 
C     A=NOISE  AMPLITUDE 
C     S=NOISE  VARIANCE 
C     AM=NOISE  MEAN 
C     AA=FILTER  COEFFICIENTS 
C     IX=SEED 

C     MU=SCALING  FACTOR  WITHIN  THE  LMS  ALGORITHM 
C     IM=INPUT  IMAGE 

C     W=WHITE  GAUSSIAN  NOISE  GENERATED  BY  SUBROUTINE 
C     U=SIGNAL  PLUS  NOISE(IM  +  W) 
C     Y=FILTER  OUTPUT 
C     E=ERROR(U  -  Y) 

C     FMINE,FMAXE,FMINU,FAMXU=PARAMETERS 
C        TO  BE  USED  TO  CONVERT  DECIMAL  DATA  TO  BYTE  DATA 
C 
C     ******* INITIAL  VALUES********* 

A=l. 

MU=35.E-9 

AA(0)=0. 

AA(1)=0. 

AA(2)=0. 

AA(3)=0. 

S=40. 

AM=0. 

IX=65539 

FMINE=1.E+10 

FMAXE=-1.E+10 

FMINU=1. E+10 

FMAXU=-1.E+10 
C 

C     ******0PEN  AN  IMAGE  FILE,  CONVERT  THE  BYTE  DATA  INTO  INTEGERS 
C     AND  THEN  PLACE  THESE  VALUES  IN  A  MATRIX******* 

OPEN  (UNIT=1,  NAME  ='H0US1G. DAT' ,  TYPE  ='OLD',  ACCESS  = 

*  'DIRECT',  REC0RDSIZE=32,  MAXREC=128) 
DO  100  M=0,127 

READ(1'M+1)  (B(J),J=0,127) 
DO  110  N=0.127 
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IF(B(N).LT.  0)THEN 
IM(M,N)=B(N)+256 
ELSE 

IM(M,N)=B(N) 
ENDIF 
110      CONTINUE 
100   CONTINUE 
C 

C     ****,v**ADD  WHITE  GAUSSIAN  NOISE  TO  THE  IMAGE  AND  SET  THE 
C     VALUES  OUTSIDE  THE  IMAGE  TO  ZERO******* 
DO  10  M=-l,127 
DO  20  N=- 1,127 

IF  ((M.  LT.  0).OR.  (N.  LT.  0))THEN 
W(M,N)=0. 
IM(M,N)=0. 
ELSE 

CALL  GAUSS(IX,S,AM,V) 
V(M,N)=V*A 
ENDIF 

U(M,N)=IM(M,N)+W(M,N) 
20       CONTINUE 
10    CONTINUE 
C 

C     ****,v*uSE  THE  LMS  ADAPTIVE  ALGORITHM  TO  UPDATE 
C     THE  FILTER  COEFICIENTS********** 
DO  40  K=0,127 
DO  50  J=0,127 

Y(K,J)=AA(0)*W(K,J)+AA(1)*W(K-1,J)+ 

*  AA(2)*W(K,J-1)+AA(3)*W(K-1,J-1) 
E(K,J)=U(K,J)-Y(K,J) 
AA(0)=AA(0)+MU*E(K,J)*W(K,J) 
AA(1)=AA( 1)+MU*E(K, J)*W(K-1, J) 
AA(2)=AA(2)+MU*E(K,J)*W(K,J-1) 
AA(3)=AA(3)+MU*E(K,J)*W(K-1,J-1) 

50       CONTINUE 
40    CONTINUE 
C 

C     ****Vfyc:SfV-cHANGE  SIGNAL  PLUS  NOISE  AND  ERROR 
C         OUTPUT  INTO  BYTE  DATA  AND  THEN  WRITE 
TO  A  FILE******** 
OPEN  ( UNIT=2 , NAME= ' ERROR.  DAT ' , TYPE= ' NEW ' , ACCESS= 

*  'DIRECT' ,REC0RDSIZE=32,MAXREC=128) 

OPEN  (UNIT=3,NAME=' SIGN0ISE.DAT' ,TYPE='NEW' ,ACCESS= 

*  'DIRECT' ,REC0RDSIZE=32,MAXREC=128) 
DO  500  1=0,127 

DO  550  J=0,127 

IF(E( I , J) .  LT.  FMINE)THEN 

FMINE=E(I,J) 
ENDIF 
IF(E(I,J).GT.  FMAXE)THEN 

FMAXE=E(I,J) 
ENDIF 
550      CONTINUE 
500   CONTINUE 

IF  (FMINE.LT.  0)  THEN 
FMAXE=FMAXE-FMINE 
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ENDIF 

DO  600  1=0,127 
DO  650  J=0,127 

IF(FMINE.LT.  0)THEN 
E(I,J)=E(I,J)-FMINE 
E(I,J)=E(I,J)''^255.  /FMAXE 
ELSE 

E( I , J)=(E( I , J)=FMINE)*255.  /FMAXE 
ENDIF 

IE=NINT(E(I,J)) 
IFdE.GT.  127)  THEN 

IE=IE-256 
ENDIF 

BYTEE(J)-IE 
650      CONTINUE 

WRITE  (2'I+1)  (BYTEE(N),N=0,127) 
600   CONTINUE 

DO  700  1=0,127 
DO  750  J=0,127 

IF(U(I,J).LT.  FMINU)THEN 

FMINU=U(I,J) 
ENDIF 
IF(U( I , J).  GT.  FMAXU)THEN 

FMAXU=U(I,J) 
ENDIF 
750      CONTINUE 
700   CONTINUE 

IF  (FMINU.  LT. 0)  THEN 

FMAXU=FMAXU-FMINU 
ENDIF 

DO  400  1=0,127 
DO  450  J=0,127 

IF(FMINU.  LT. 0)THEN 
U(I,J)=U(I,J)-FMINU 
U(  I ,  J)=U(  I ,  J)'V255.  /FMAXU 
ELSE 

U(  I ,  J)=(U(  I ,  J)=FMINU)'V255.  /FMAXU 
ENDIF 

IU=NINT(U(I,J)) 
IF(IU.  GT.  127)  THEN 

IU=IU-256 
ENDIF 

BYTEU(J)-IU 
450      CONTINUE 

WRITE  (3'I+1)  (BYTEU(N),N=0,127) 
400   CONTINUE 

CLOSE  (UNIT=2) 
CLOSE  (UNIT=3) 
STOP 
END 
C 

SUBROUTINE  GAUSS(IX,S,AM,V) 

AMP=0. 0 

DO  50  1=1,12 

RANDOM=RAN( IX) 
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AMP=AMP+RANDOM 
50    CONTINUE 

V=(AMP-6.0)*S+AM 

RETURN 

END 
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APPENDIX  C.     LMS  ALGORITHM  IMPLEMENTED  IN  AN  ADAPTIVE 

LINE  ENHANCER 

C  THIS   IS  A  VAX/ VMS  FORTRAN  PROGRAM  THAT  IMPLEMENTS  A  TWO 

C  DIMENSIONAL  2X2  ADAPTIVE  LMS  FILTER  WITHIN  AN  ADAPTIVE 

C  LINE  ENHANCER 

C 

INTEGER  M,N,K,J 

BYTE   B(0:  127) ,BYTEY(0:  127) ,BYTEU(0: 127) 

INTEGER*4   INTY(0:  127,0: 127),INTU(0:  127,0:  127),IY,IU 

REAL*4  MU,A,FMINY,FMAXY,FMINU,FMAXU 

REAL*4  AA(0:3>,E(0:  127,0:  127), U(0:  127,0:  127) 

REAL*4  W(-l:  127,-1:  127), Y(0:  127,0:  127),IM(0:  127,0:  127) 
C 

C     *,v,v,v,v,vvARIABLE  DEFINITIONS****** 
C     A=NOISE  AMPLITUDE 
C     S=NOISE  VARIANCE 
C     AM=NOISE  MEAN 
C     AA=FILTER  COEFFICIENTS 
C     IX=SEED 

C     MU=SCALING  FACTOR  WITHIN  THE  LMS  ALGORITHM 
C     IM=INPUT  IMAGE 

C     W=WHITE  GAUSSIAN  NOISE  GENERATED  BY  SUBROUTINE 
C     U=SIGNAL  PLUS  NOISE(IM  +  W) 
C     Y=FILTER  OUTPUT 

C     WW=DELAYED  VERSION  OF  SIGNAL  PLUS  NOISE(U) 
C     E=ERROR(U  -  Y) 

C     FMINY,FMAXY,FMINU,FAMXU=PARAMETERS 
C        TO  BE  USED  TO  CONVERT  DECIMAL  DATA  TO  BYTE  DATA 
C 
C     ******* INITIAL  VALUES********* 

A=l. 

MU=35.E-9 

AA(0)=0. 

AA(1)=0. 

AA(2)=0. 

AA(3)=0. 

S=40. 

AM=0. 

IX=65539 

FMINU=1.  E+10 

FMAXU=-1.E+10 

FMINY=1.  E+10 

FMAXY=-1.E+10 
C 

C     ******OPEN  AN  IMAGE  FILE,  CONVERT  THE  BYTE  DATA  INTO  INTEGERS 
C     AND  THEN  PLACE  THESE  VALUES  IN  A  MATRIX******* 

OPEN  (UNIT=1,  NAME  ='H0US1G. DAT' ,  TYPE  ='OLD' ,  ACCESS  = 
*       'DIRECT',  REC0RDSIZE=32,  MAXREC=128) 

DO  100  M=0,127 

READ(1'M+1)  (B(J),J=0,127) 
DO  110  N=0,127 
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IF(B(N).LT.  0)THEN 
IM(M,N)=B(N)+256 
ELSE 

IM(M,N)=B(N) 
END  IF 
110      CONTINUE 
100   CONTINUE 
C 

C     ******* ADD  WHITE  GAUSSIAN  NOISE  TO  THE  IMAGE  AND  SET  THE 
C     VALUES  OUTSIDE  THE  IMAGE  TO  ZERO******* 
DO  10  M=- 1,127 
DO  20  N=-l,127 

IF  ((M.  LT.  0).OR.  (N.  LT.  0))THEN 
W(M,N)=0. 
IM(M,N)=0. 
ELSE 

CALL  GAUSS(IX,S,AM,V) 
W(M,N)=V*A 
ENDIF 

U(M,N)=IM(M,N)+W(M,N) 
WW(M,N)=U(M-1,N) 
20        CONTINUE 
10    CONTINUE 
C 

C     ******USE  THE  LMS  ADAPTIVE  ALGORITHM  TO  UPDATE 
C     THE  FILTER  COEFICIENTS********** 
DO  40  K=0,127 
DO  50  J=0,127 

Y(K,J)=AA(0)*WV(K,J)+AA(1)*WW(K-1,J)+ 

*  AA(2)*WW(K,J-1)+AA(3)*WW(K-1,J-1) 
E(K,J)=U(K,J)-Y(K,J) 
AA(0)=AA(0)+MU*E(K,J)*WW(K,J) 

AA( 1)=AA( 1)+MU*E(K, J)*WW(K-1, J) 

AA(2)=AA(2)+MU*E(K,J)*WV(K,J-1) 

AA(3)=AA(3)+MU*ECK,J)*WV(K-1,J-1) 
50       CONTINUE 
40    CONTINUE 
C 

C     ********CHANGE  SIGNAL  PLUS  NOISE  AND  FILTER  OUTPUT 
C         INTO  BYTE  DATA  AND  THEN  WRITE  TO  A  FILE****** 

OPEN  (UNIT=3,NAME=' SIGNOISE.DAT' ,TYPE='NEW' ,ACCESS= 

*  'DIRECT' ,REC0RDSIZE=32,MAXREC=128) 

OPEN  (UNIT=4,NAME=' FILTERED. DAT' ,TYPE='NEW' ,ACCESS= 

*  'DIRECT' ,REC0RDSIZE=32,MAXREC=128) 
DO  700  1=0,127 

DO  750  J=0,127 

IF(U(I,J).LT.  FMINU)THEN 

FMINU=U(I,J) 
ENDIF 
IF(U( I , J). GT. FMAXU)THEN 

FMAXU=U(I,J) 
ENDIF 
750      CONTINUE 
700   CONTINUE 

IF  (FMINU. LT.  0)  THEN 
FMAXU=FMAXU-FMINU 
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END  IF 

DO  400  1=0,127 
DO  450  J=0,127 

IFCFMINU. LT.  0)THEN 
U(I,J)=U(I,J)-FMINU 
U( I , J)=U( I , J)*255.  /FMAXU 
ELSE 

U( I , J)=(U( I , J)=FMINU)*255. /FMAXU 
ENDIF 

IU=NINT(U(I,J)) 
IF(IU.  GT.  127)  THEN 

IU=IU-256 
ENDIF 

BYTEU(J)-IU 
450      CONTINUE 

WRITE  (3'I+1)  (BYTEU(N),N=0,127) 
400   CONTINUE 

DO  800  1=0,127 
DO  850  J=0,127 

IF( Y( I , J).  LT.  FMINY)THEN 

FMINY=Y(I,J) 
ENDIF 
IF(Y(I,J).GT.  FMAXY)THEN 

FMAXY=Y(I,J)         — 
ENDIF 
850      CONTINUE 
800   CONTINUE 

IF  (FMINY. LT. 0)  THEN 

FMAXY=FMAXY-FMINY 
ENDIF 

DO  900  1=0,127 
DO  950  J=0,127 

IF(FMINY.  LT.  0)THEN 
Y(I,J)=Y(I,J)-FMINY 
Y( I , J)=Y( I , J)*255.  /FMAXY 
ELSE 

Y( I , J)=(Y( I , J)=FMINY)*255.  /FMAXY 
ENDIF 

IY=NINT(Y(I,J)) 
IF(IY. GT. 127)  THEN 

IY=IY-256 
ENDIF 

BYTEY(J)-IY 
950      CONTINUE 

WRITE  (4'I+1)  (BYTEY(N),N=0,127) 
900   CONTINUE 

CLOSE  (UNIT=3) 
CLOSE  (UNIT=4) 
STOP 
END 
C 

SUBROUTINE  GAUSS( IX,S,AM,V) 

AMP=0.  0 

DO  50   1=1,12 

RANDOM=RAN( IX) 
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AMP=AMP+RANDOM 
50    CONTINUE 

V=(AMP-6.0)'''S+AM 

RETURN 

END 
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APPENDIX  D.     RLS  SYSTEM  IDENTIFICATION 

C     SYSTEM  IDENTIFICATION  OF  AN  ADAPTIVE  RLS  2X2  FILTER  VERSUS 
C         2X2  RLS  FILTER  VERSUS  A  A  2X2  FILTER  WITH 
C         KNOWN  COEFFICIENTS 
C 

DIMENSION  W(-l:  32,-1:  32),P(0:  3,0:  3),AA(0:  3),X(0:  3) 
C 

C     ******VARIABLE  DEFINITIONS****** 
C     A=NOISE  AMPLITUDE 
C     S=NOISE  VARIANCE 
C     AM=NOISE  MEAN 
C     IX=SEED 

C     AA=FILTER  COEFFICIENTS 

C     W=WHITE  GAUSSIAN  NOISE  GENERATED  BY  SUBROUTINE 
C     D=DESIRED  OUTPUT 
C     Y=FILTER  OUTPUT 
C     E=ERROR(D  -  Y) 
C     MXM  =  ARRAY  SIZE 
C     NXN  =  FILTER  SIZE 
C     Q=WEIGHTING  FACTOR 
C     P=INVERSE  OF  COVARIANCE  MATRIX 
C 
C     -itif/Hticit*   INPUT  INITIAL  VALUES  ******i(***** 

M=32 

N=2 

A=l. 

AA(0)=0. 

AA(1)=0. 

AA(2)=0. 

AA(3)=0. 

Q=l. 
S=l 
AM=0 
IX=65539 
C 

C     Vr-vV"V:t*  BUILD  INITIAL  'P'  MATRIX  ********* 
DO  1  K=0,(N**2-1) 
DO  5  J=0,(N**2-1) 
IF(K.  EQ.  J)THEN 
P(K,J)=100. 
ELSE 

P(K,J)=0.0 
END  IF 
5        CONTINUE 
1     CONTINUE 

OPEN  (UNIT=63,FILE='RLSID2D2  DATA ', STATUS= ' OLD ' ) 
C 

C     *******  CALCULATE  INPUT  MATRIX  (WHITE  GAUSSIAN  NOISE)  ****** 
DO  10  K=-1,M-1 
DO  20  J=-1,M-1 

IF((K.  LT.  0).OR.  (J.  LT.  0))THEN 
W(K,J)=0. 
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ELSE 

CALL  GAUSS(IX,S,AM,V) 
W(K,J)=V'^A 
END  IF 
20       CONTINUE 
10    CONTINUE 
C 

C     **,v****Vf,v  CALCULATE  ERROR  BETWEEN  ADAPTIVE  FILTER  OUTPUT 
C         AND  KNOWN  FILTER  OUTPUT.   USE  THIS  ERROR  TO  UPDATE 
C         FILTER  COEFFICIENTS  AND  THE  'P'  MATRIX  ******^v*yr**** 
DO  40  K=0,M-1 
DO  50  J=0,M-1 
L=K*M+J 

D=. 4*W(K,J)+. 6*W(K-1,J)-.  3*W(K,J-1)+.  2*W(K-1,J-1) 
Y=AA(0)*W(K,J)+AA(1)*W(K-1,J)+AA(2)*W(K,J-1)+ 

*  AA(3)*W(K-1,J-1) 
E=D-Y 

X(0)=W(K,J) 
X(1)=W(K-1,J) 
X(2)=W(K,J-1) 
X(3)=W(K-1,J-1) 

CALL  RLS(P,X,AA,N,E,Q) 
PRINT  15,  L,E,AAC0),AA(1),AA(2),AA(3) 
WRITE  (63,15)  L,E,AA(0),AA(1),AA(2),AA(3) 
50       CONTINUE 
40    CONTINUE 

15    FORMAT  ('  ' ,2X,I3,2X,F9.5,2X,F9.5,2X,F9.5,2X,F9.5,2X,F9.5) 
STOP 
END 
C 

c 

C     V!"jV*v«v**'jV'>v  to  COMPUTE  THE  GAIN  MATRIX  ********** 
SUBROUTINE  RLS(P,X,AA,N,E,Q) 
DIMENSION  P(0:  3,0:  3)  ,X(0:  3),AA(0:  3),TEMY(0:  3,0:  3), 

*  GAMA(0:3),TEMC(0:  3) 
DO  30  L=0,3 

TEM=0. 
DO  20  K=0,3 
20       TEM=TEM+X(K)*P(K,L) 
30    GAMA(L)=TEM 
GAM=Q 

DO  40  K=0,3 
40    GAM=GAM+GAMA(K)*X(K) 
DO  60  L=0,3 
TEM=0. 

DO  50  K=0,3 
50       TEM=TEM+P(L,K)*X(K) 
60    TEMC(L)=TEM 
DO  70  L=0,3 
DO  70  K=0,3 
70    TEMY(L,K)=TEMC(L)*GAMA(K) 
DO  80  K=0,3 
DO  80  L=0,3 
80    P(K,L)=(P(K,L)-(TEMY(K,L)/GAM))/Q 
DO  100  K=0,3 
TEM=0. 
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DO  90  L=0,3 
90       TEM=TEM+P(K,L)*X(L) 
100   TEMC(K)=TEM 

DO  110  K=0,3 
110   AA(K)=AA(K)+TEMC(K)*E 

Q=.  99''^Q+.  01 

RETURN 

END 
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APPENDIX  E.     RLS  ALGORITHM  IMPLEMENTED  IN  A  NOISE 

CANCELER 

C     THIS  IS  A  VAX/VMS  FORTRAN  PROGRAM  THAT  IMPLEMENTS  A  TWO 
C     DIMENSIONAL  2X2  ADAPTIVE  RLS  FILTER  WITHIN  A  NOISE 
C     CANCELER 
C 

INTEGER  M,N,K,J 

BYTE  B(0:  127) ,BYTEE(0: 127) ,BYTEU(0: 127) 

INTEGER*4  INTE(0:  127,0:  127),INTU(0:  127,0:  127),IE,IU 

REAL*4  A,FMINE,FMAXE,FMINU,FMAXU 

REAL*4  AA(0: 3),E(0: 127,0:  127), U(0:  127,0:  127) 

REAL*4  P(0:  3,0:  3),X(0:  3),IM(0:  127,0:  127) 

REAL*4  W(-l:  127,-1:  127), Y(0:  127,0:  127) 
C 

C     ******VARIABLE  DEFINITIONS****** 
C     A=NOISE  AMPLITUDE 
C     S=NOISE  VARIANCE 
C     AM=NOISE  MEAN 
C     AA=FILTER  COEFFICIENTS 
C     IX=SEED 
C     I M= INPUT  IMAGE 

C     W=WHITE  GAUSSIAN  NOISE  GENERATED  BY  SUBROUTINE 
C     U=SIGNAL  PLUS  NOISE(IM  +  W) 
C     Y=FILTER  OUTPUT 
C     E=ERROR(U  -  Y) 
C     Q=WEIGHTING  FACTOR 
C     P=INVERSE  OF  COVARIANCE  MATRIX 
C     FMINE , FMAXE , FMINU , FAMXU=PARAMETERS 
C        TO  BE  USED  TO  CONVERT  DECIMAL  DATA  TO  BYTE  DATA 
C     MXM=ARRAY  SIZE 
C     NXN=FILTER  SIZE 
C 
C     ******* INITIAL  VALUES********* 

M=128 

N=2 

Q=l. 

A=l. 

AA(0)=0. 

AA(1)=0. 

AA(2)=0. 

AA(3)=0. 

S=40. 

AM=0. 

IX=65539 

FMINE=1. E+10 

FMAXE=-1.E+10 

FMINU=1.E+10 

FMAXU=-1.E+10 
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c 

C     ******0PEN  AN  IMAGE  FILE,  CONVERT  THE  BYTE  DATA  INTO  INTEGERS 
C     AND  THEN  PLACE  THESE  VALUES  IN  A  MATRIX**'^*'^?* 

OPEN  (UNIT=1,  NAME  =' HOUSIG. DAT* ,  TYPE  ='OLD',  ACCESS  = 

*  'DIRECT',  REC0RDSIZE=32,  MAXREC=128) 
DO  100  K=0,127 

READ(1'K+1)  (B(L),L=0,127) 
DO  110  J=0,127 

IF(B(J).LT.  0)THEN 
IM(K,J)=B(J)+256 
ELSE 

IM(K,J)=B(J) 
END  IF 
110      CONTINUE 
100   CONTINUE 

CLOSE  (UNIT=1) 
C 

C     ******  BUILD  INITIAL  'P'  MATRIX  ********* 
DO  1  K=0,(N**2-1) 
DO  5  J=0,(N**2-1) 
IF(K.  EQ.  J)THEN 
P(K,J)=100. 
ELSE 

P(K,J)=0.  0 
ENDIF  ^ 

5        CONTINUE 
1     CONTINUE 
C 

C     *******ADD  WHITE  GAUSSIAN  NOISE  TO  THE  IMAGE  AND  SET  THE 
C     VALUES  OUTSIDE  THE  IMAGE  TO  ZERO******* 
DO  10  K=-1,M-1 
DO  20  J=-1,M-1 

IF  ((K.  LT.  0).OR.  (J.  LT.  0))THEN 
W(K,J)=0. 
IM(K,J)=0. 
ELSE 

CALL  GAUSS(IX,S,AM,V) 
V(K,J)=V'-A 
ENDIF 

U(K,J)=IM(K,J)+V(K,J) 
20        CONTINUE 
10    CONTINUE 
C 

C     *********  CALCULATE  ERROR  BETWEEN  ADAPTIVE  FILTER  OUTPUT 
C         AND  KNOWN  FILTER  OUTPUT.   USE  THIS  ERROR  TO  UPDATE 
C         FILTER  COEFFICIENTS  AND  THE  'P'  MATRIX  ************* 
DO  40  K=0,M-1 
DO  50  J=0,M-1 
L=K*M+J 
Y=AA(0)*W(K,J)+AA(1)*W(K-1,J)+AA(2)*W(K,J-1)+ 

*  AA(3)*W(K-1,J-1) 
E(K,J)=U(K,J)-Y(K,J) 
EE=E(K,J) 
X(0)=W(K,J) 
X(1)=W(K-1,J) 
X(2)=W(K,J-1) 
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X(3)=W(K-1,J-1) 

CALL  RLS(P,X,AA,N,EE,Q) 
50       CONTINUE 
40    CONTINUE 
C 

C     ***,v*,v,v,vcHANGE  SIGNAL  PLUS  NOISE  AND  ERROR  OUTPUT 
C         INTO  BYTE  DATA  AND  WRITE  TO  A  FILE***** 

OPEN  (UNIT=2,NAME=' ERROR.  DAT' ,TYPE='NEW' ,ACCESS= 

*  'DIRECT '.RECORDS I ZE=32,MAXREC=128) 

OPEN  (UNIT=3.NAME=' SIGNOISE.DAT' ,TYPE='NEW' ,ACCESS= 

*  •DIRECT\REC0RDSIZE=32,MAXREC=128) 
DO  500   1=0,127 

DO  550  J=0,127 

IF(E(I,J).LT.  FMINE)THEN 

FMINE=E(I,J) 
END  IF 
IF(E( I , J). GT.  FMAXE)THEN 

FMAXE=E(I,J) 
ENDIF 
550      CONTINUE 
500   CONTINUE 

IF  (FMINE. LT.  0)  THEN 

FMAXE=FMAXE-FMINE 
ENDIF 

DO  600  1=0,127 
DO  650  J=0,127 

IF(FMINE.  LT. 0)THEN 
E(I,J)=E(I,J)-FMINE 
E( I , J)=E(I , J)*255. /FMAXE 
ELSE 

E( I , J)=(E( I , J)=FMINE)*255. /FMAXE 
ENDIF 

IE=NINT(E(I,J)) 
IFdE.GT.  127)  THEN 

IE=IE-256 
ENDIF 

BYTEE(J)-IE 
650      CONTINUE 

WRITE  (2'I+1)  (BYTEE(N),N=0,127) 
600   CONTINUE 

DO  700  1=0,127 
DO  750  J=0,127 

IF(U(I,J).LT.  FMINU)THEN 

FMINU=U(I,J) 
ENDIF 
IF(U( I , J).  GT. FMAXU)THEN 

FMAXU=U(I,J) 
ENDIF 
750      CONTINUE 
700   CONTINUE 

IF  (FMINU. LT. 0)  THEN 

FMAXU=FMAXU-FMINU 
ENDIF 

DO  400   1=0,127 
DO  450  J=0,127 

IF(FMINU. LT. 0)THEN 
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U(I,J)=U(I,J)-FMINU 
U(I,J)=U(I,J)*255.  /FMAXU 
ELSE 

U( I , J)=(U( I , J)=FMINU)*255.  /FMAXU 
END  IF 

IU=NINT(U(I,J)) 
IF(IU.  GT.  127)  THEN 

IU=IU-256 
ENDIF 

BYTEU(J)-IU 
450      CONTINUE 

WRITE  (3'I+1)  (BYTEU(N),N=0,127) 
400   CONTINUE 

CLOSE  (UNIT=2) 
CLOSE  (UNIT=3) 
STOP 
END 
C 

C     **'sV')V*****  jQ   COMPUTE  THE  GAIN  MATRIX  ********** 
SUBROUTINE  RLS(P,X, AA,N,E ,Q) 

DIMENSION  P(0:  3,0:  3),X(0:  3),AA(0:  3),TEMY(0:  3,0:  3), 
*      GAMA(0:  3),TEMC(0:  3) 
DO  30  L=0,3 
TEM=0. 

DO  20  K=0,3 
20       TEM=TEM+X(K)*P(K,L) 
30    GAMA(L)=TEM 
GAM=Q 

DO  40  K=0,3 
40    GAM=GAM+GAMA(K)*X(K) 
DO  60  L=0,3 
TEM=0. 

DO  50  K=0,3 
50       TEM=TEM+P(L,K)*X(K) 
60    TEMC(L)=TEM 
DO  70  L=0,3 
DO  70  K=0,3 
70    TEMY(L,K)=TEMC(L)*GAMA(K) 
DO  80  K=0,3 
DO  80  L=0,3 
80    P(K,L)=(P(K,L>-(TEMY(K,L)/GAM))/Q 
DO  100  K=0,3 
TEM=0. 
DO  90  L=0,3 
90       TEM=TEM+P(K,L)*X(L) 
100   TEMC(K)=TEM 

DO  110  K=0,3 
110   AA(K)=AA(K)+TEMC(K)*E 
Q=.  99*Q+.  01 
RETURN 
END 


58 


APPENDIX  F.     RLS  ALGORITHM  IMPLEMENTED  IN  AN  ADAPTIVE 


C 

c 
c 
c 


c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


LINE  ENHANCER 

THIS  IS  A  VAX/ VMS  FORTRAN  PROGRAM  THAT  IMPLEMENTS  A  TWO 
DIMENSIONAL  2X2  ADAPTIVE  RLS  FILTER  WITHIN  AN  ADAPTIVE 
LINE  ENHANCER 

INTEGER  M,N,K,J 

BYTE  B(0:  127) ,BYTEU(0: 127) ,BYTEY(0: 127) 

INTEGER*4  INTY(0:  127,0:  127),INTU(0:  127,0:  127),IU,IY 

REAL*4  A,FMINU,FMAXU,FMINY,FMAXY 

REAL*4  AA(0: 3),E(0: 127,0: 127), U(0: 127,0:  127) 

REAL*4  P(0:3,0:3),X(0:3),IM(0:  127,0:  127) 

REAL*4  W(-l: 127,-1: 127), Y(0: 127,0: 127) 

*****,vvARIABLE  DEFINITIONS****** 

A=NOISE  AMPLITUDE 

S=NOISE  VARIANCE 

AM=NOISE  MEAN 

AA=FILTER  COEFFICIENTS 

IX=SEED 

I M= INPUT  IMAGE 

W=WHITE  GAUSSIAN  NOISE  GENERATED  BY  SUBROUTINE 

U=SIGNAL  PLUS  NOISE(IM  +  W) 

Y=FILTER  OUTPUT 

E=ERROR(U  -  Y) 

WV=DELAYED  VERSION  OF  SIGNAL  PLUS  NOISE(U) 

Q=WEIGHTING  FACTOR 

P=INVARSE  OF  COVARIANCE  MATRIX 

FMINU , FAMXU , FMINY , FMAXY=PARAMETERS 

TO  BE  USED  TO  CONVERT  DECIMAL  DATA  TO  BYTE  DATA 
MXM=ARRAY  SIZE 
NXN=FILTER  SIZE 


m3m  Jm  KM  9j0*J^  mi 


'-** INITIAL  VALUES********* 
M=128 
N=2 
Q=l. 
A=l. 

AA(0)=0. 
AA(1)=0. 
AA(2)=0. 
AA(3)=0. 
S=40. 
AM=0. 
IX=65539 
FMINU=1.E+10 
FMAXU=-1.E+10 
FMINY=1.E+10 
FMAXY=-1.E+10 
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c 
c 


110 
100 

c 
c 


5 
1 
C 
C 

c 


20 
10 
C 

c 
c 


46 
45 
C 
C 


***,v,vyroPEN  AN  IMAGE  FILE,  CONVERT  THE  BYTE  DATA  INTO  INTEGERS 
AND  THEN  PLACE  THESE  VALUES  IN  A  MATRIX^"****** 
OPEN  (UNIT=1,  NAME  ='H0US1G. DAT' ,  TYPE  ='OLD',  ACCESS  = 
*       'DIRECT',  RECORDS I ZE=32,  MAXREC=128) 
DO  100  K=0,127 

READ(1'K+1)  (B(L),L=0,127) 
DO  110  J=0,127 

IF(B(J).LT.  0)THEN 
IM(K,J)=B(J)+256 
ELSE 

IM(K,J)=B(J) 
ENDIF 
CONTINUE 
CONTINUE 
CLOSE  (UNIT=1) 

******  BUILD  INITIAL  'P'  MATRIX  ********* 
DO  1  K=0,(N**2-1) 
DO  5  J=0,(N**2-1) 
IF(K.  EQ.  J)THEN 
P(K,J)=100. 
ELSE 

P(K,J)=0.  0 
ENDIF 
CONTINUE 
CONTINUE 

*******ADD  WHITE  GAUSSIAN  NOISE  TO  THE  IMAGE  AND  SET  THE 
VALUES  OUTSIDE  THE  IMAGE  TO  ZERO******* 
DO  10  K=-2,M-1 
DO  20  J=-2,M-1 

IF  ((K.  LT.  0).OR.  (J.  LT.  0))THEN 
W(K,J)=0. 
IM(K,J)=0. 
ELSE 

CALL  GAUSS(IX,S,AM,V) 
W(K,J)=V*A 
ENDIF 

U(K,J)=IM(K,J)+W(K,J) 
CONTINUE 
CONTINUE 

*****COMPUTE  THE  DELAYED  VERSION  OF  THE  SIGNAL 

PLUS  NOISE  MATRIX(U)****** 
DO  45  K=-l,127 
DO  46  J=- 1,127 

IF((M.  LT. -1).0R.  (N.  LT. -1)  THEN 

WW(K,J)=0. 
ELSE 

WW(K,J)=U(K-1,J) 
ENDIF 
CONTINUE 
CONTINUE 

*********  CALCULATE  ERROR  BETWEEN  ADAPTIVE  FILTER  OUTPUT 
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50 

40 

C 

C 

C 


C         AND  KNOWN  FILTER  OUTPUT.   USE  THIS  ERROR  TO  UPDATE 
C         FILTER  COEFFICIENTS  AND  THE  'P'  MATRIX  *Vr**yc*Tr,v*,v*** 
DO  40  K=0,M-1 
DO  50  J=0,M-1 

Y=AA(0)'VWV(K,J)+AA(1)*WV(K-1,J)+AA(2)*WW(K,J-1)+ 

*  AA(3)*WW(K-1,J-1) 
E(K,J)=U(K,J)-Y(K,J) 
EE=E(K,J) 
X(0)=WW(K,J) 
X(1)=WW(K-1,J) 
X(2)=WV(K,J-1) 
X(3)=WV(K-1,J-1) 
CALL  RLS(P,X,AA,N,EE,Q) 

CONTINUE 
CONTINUE 

********CHANGE  SIGNAL  PLUS  NOISE  AND  FILTER  OUTPUT 

INTO  BYTE  DATA  AND  WRITE  TO  A  FILE***** 
OPEN  (UNIT=3.NAME=' SIGNOISE.DAT' ,TYPE='NEW' ,ACCESS= 

*  'DIRECT^ ,REC0RDSIZE=32,MAXREC=128) 

OPEN  (UNIT=4NAME=' FILTERED.  DAT' ,TYPE='NEW' ,ACCESS= 

*  'DIRECT' ,REC0RDSIZE=32,MAXREC=128) 
DO  700  1=0,127 

DO  750  J=0,127 

IF(U(I,J).LT.  FMINU)THEN 

FMINU=U(I,J) 
ENDIF 
IF( U( I , J) . GT. FMAXU)THEN 

FMAXU=U(I,J) 
ENDIF 
750      CONTINUE 
700   CONTINUE 

IF  (FMINU.  LT.  0)  THEN 

FMAXU=FMAXU-FMINU 
ENDIF 

DO  400  1=0,127 
DO  450  J=0,127 

IF(  FMINU.  LT.O) THEN 
U(I,J)=U(I,J)-FMINU 
U( I , J)=U( I , J)*255.  /FMAXU 
ELSE 

U(I,J)=(U(I,J)=FMINU)*255. /FMAXU 
ENDIF 

IU=NINT(U(I,J)) 
IF(IU.  GT.  127)  THEN 

IU=IU-256 
ENDIF 

BYTEU(J)-IU 
450      CONTINUE 

WRITE  (3'I+1)  (BYTEU(N),N=0,127) 
400   CONTINUE 

DO  800  1=0,127 
DO  850  J=0,127 

IF(Y(I,J).LT.  FMINY)THEN 

FMINY=Y(I,J) 
ENDIF 
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IF( Y( I , J). GT. FMAXY)THEN 

FMAXY=Y(I,J) 
ENDIF 
850      CONTINUE 
800   CONTINUE 

IF  (FMINY. LT. 0)  THEN 

FMAXY=FMAXY-FMINY 
ENDIF 

DO  900  1=0,127 
DO  950  J=0,127 

IF(FMINY. LT. 0)THEN 
Y(I,J)=Y(I,J)-FMINY 
Y(I,J)=Y(I,J)*255.  /FMAXY 
ELSE 

Y(I,J)=(Y(I,J)=FMINY)*255.  /FMAXY 
ENDIF 

IY=NINT(Y(I,J)) 
IFdY.GT.  127)  THEN 

IY=IY-256 
ENDIF 

BYTEY(J)-IY 
950      CONTINUE 

WRITE  (4'I+1)  (BYTEY(N),N=0,127) 
900   CONTINUE 

CLOSE  (UNIT=3) 
CLOSE  (UNIT=4) 
STOP 
END 

*****Vf**vr  TO  COMPUTE  THE  GAIN  MATRIX  ********Vf* 
SUBROUTINE  RLS(P,X,AA,N,E,Q) 

DIMENSION  P(0:  3,0:  3)  ,X(0:  3)  ,AA(0:  3)  ,TEMY(0:  3,0:  3), 
*      GAMA(0:3),TEMC(0:3) 
DO  30  L=0,3 
TEM=0. 
DO  20  K=0,3 
20       TEM=TEM+X(K)^^P(K,L) 
30    GAMA(L)=TEM 
GAM=Q 

DO  40  K=0,3 
40    GAM=GAM+GAMA(K)*X(K) 
DO  60  L=0,3 
TEM=0. 

DO  50  K=0,3 
50       TEM=TEM+P(L,K)*X(K) 
60    TEMC(L)=TEM 
DO  70  L=0,3 
DO  70  K=0,3 
70    TEMY(L,K)=TEMC(L)*GAMA(K) 
DO  80  K=0,3 
DO  80  L=0,3 
80    P(K,L)=(P(K,L)-(TEMY(K,L)/GAM))/Q 
DO  100  K=0,3 
TEM=0. 


C 
C 
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DO  90  L=0,3  APP02140 

90       TEM=TEM+P(K,L)*X(L)  APP02150 

100   TEMC(K)=TEM  APP02160 

DO  110  K=0,3  APP02170 

110   AA(K)=AA(K)+TEMC(K)*E  APP02180 

Q=.  99'VQ+.  01  APP02190 

RETURN  APP02200 

END  APP02210 
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